
# Replication file for Appendix G.2., Table 17, 19, 21, 23, 25


#install.packages("Hmisc")
require(Hmisc)
#install.packages("stargazer")
require(stargazer)
require(multiwayvcov)
require(lmtest)

load("Brexit_Rep.RData")

treated.sub <-  brexit.sub2[which(brexit.sub2$post == 1),]
control.sub <-  brexit.sub2[which(brexit.sub2$post == 0),]


varlabels <- c("Treat","Constant")

right.side.eq <- "~ wave"

### Control Wave 6 and Control Wave 7 (June 2016, Pre Referendum/ Nov. 2016)

ones <- rep_len(1, length.out = nrow(control.sub))
twos <- rep_len(2, length.out = nrow(control.sub))


dvs.6 <- cbind(ones, control.sub$w8w6, control.sub$leaver, control.sub$migs.take.jobs.6, control.sub$migs.more.terror.6,control.sub$closed.immigrants.6,
               control.sub$hurt.standing.refugee.6,
               control.sub$threaten.culture.refugee.6,control.sub$overwhelm.welfare.refugee.6, control.sub$GUNQID)

dvs.7 <- cbind(twos,control.sub$w8w6, control.sub$leaver, control.sub$migs.take.jobs.7,control.sub$migs.more.terror.7,control.sub$closed.immigrants.7,
               control.sub$hurt.standing.refugee.7,
               control.sub$threaten.culture.refugee.7,control.sub$overwhelm.welfare.refugee.7, control.sub$GUNQID)

dvs <- Cs(migs.take.jobs,migs.more.terror,closed.immigrants,hurt.standing.refugee,
                threaten.culture.refugee,overwhelm.welfare.refugee)


c67 <- as.data.frame(rbind(dvs.6, dvs.7))


names(c67) <- c("wave", "w8w6", "leaver",dvs,"GUNQID")

c67.l <- subset(c67, leaver == 1)
c67.r <- subset(c67, leaver == 0)

lf.l <- vector(length(dvs), mode = "list")
names(lf.l) <- dvs

lf.r <- vector(length(dvs), mode = "list")
names(lf.r) <- dvs


for (i in 1:(length(dvs))){
  
  modelformula <- paste(dvs[i],right.side.eq)
  print(modelformula)
  
  x <- eval(substitute(lm(.modelformula, data = c67.l, 
                                        weights = w8w6), list(.modelformula = modelformula)))
  

  
 x.clust <- cluster.vcov(x, c67.l$GUNQID)
 x.mcse <- coeftest(x, x.clust)
  
  
lf.l[[dvs[i]]] <- x.mcse
  
print(lf.l[[dvs[i]]])

}


for (i in 1:(length(dvs))){
  
  modelformula <- paste(dvs[i],right.side.eq)
  print(modelformula)
  
  x.r <- eval(substitute(lm(.modelformula, data = c67.r, 
                                       weights = w8w6), list(.modelformula = modelformula)))
  
  x.clust.r <- cluster.vcov(x.r, c67.r$GUNQID)
  x.mcse.r <- coeftest(x.r, x.clust.r)
  
  lf.r[[dvs[i]]] <- x.mcse.r
  
  print(lf.r[[dvs[i]]])
  
}


# Refugees

stargazer(lf.l[[dvs[5]]], lf.r[[dvs[5]]], lf.l[[dvs[6]]], lf.r[[dvs[6]]],lf.l[[dvs[4]]], lf.r[[dvs[4]]], covariate.labels = varlabels)

# Migrants

stargazer(lf.l[[dvs[3]]], lf.r[[dvs[3]]], lf.l[[dvs[1]]], lf.r[[dvs[1]]],lf.l[[dvs[2]]], lf.r[[dvs[2]]], covariate.labels = varlabels)


#### Treated Wave 5 and Wave 6  (Dec. 2015/ July 2016, Post Referendum)

ones <- rep_len(1, length.out = nrow(treated.sub))
twos <- rep_len(2, length.out = nrow(treated.sub))

t.ref.test.5 <- cbind(ones,treated.sub$w8w6,treated.sub$leaver,treated.sub$hurt.standing.refugee.5,
                      treated.sub$threaten.culture.refugee.5,treated.sub$overwhelm.welfare.refugee.5, 
                      treated.sub$closed.immigrants.5, treated.sub$GUNQID)
t.ref.test.6 <- cbind(twos,treated.sub$w8w6,treated.sub$leaver,treated.sub$hurt.standing.refugee.6,
                      treated.sub$threaten.culture.refugee.6,treated.sub$overwhelm.welfare.refugee.6,
                      treated.sub$closed.immigrants.6, treated.sub$GUNQID)

t.w56 <- as.data.frame(rbind(t.ref.test.5, t.ref.test.6))

ref.names <- c("hurt.standing.refugee","threaten.culture.refugee","overwhelm.welfare.refugee","closed.immigrants")

names(t.w56) <- c("wave", "w8w6", "leaver",ref.names, "GUNQID")

t.w56.l <- subset(t.w56, leaver == 1)
t.w56.r <- subset(t.w56, leaver == 0)

tl <- vector(length(ref.names), mode = "list")
names(tl) <- ref.names

tr <- vector(length(ref.names), mode = "list")
names(tr) <- ref.names

#Leavers
for (i in 1:(length(ref.names))){
  
  modelformula <- paste(ref.names[i],right.side.eq)
  print(modelformula)
  
  x.l <- eval(substitute(lm(.modelformula, data = t.w56.l, 
                                       weights = w8w6), list(.modelformula = modelformula)))
  
  
  
  x.clust.l <- cluster.vcov(x.l, t.w56.l$GUNQID)
  x.mcse.l <- coeftest(x.l, x.clust.l)
  
  
  tl[[ref.names[i]]] <- x.mcse.l
  
  print(tl[[ref.names[i]]])
  
}

#Remainers

for (i in 1:(length(ref.names))){
  
  modelformula <- paste(ref.names[i],right.side.eq)
  print(modelformula)
  
 x.r <- eval(substitute(lm(.modelformula, data = t.w56.r, 
                                       weights = w8w6), list(.modelformula = modelformula)))
  
  
 x.clust.r <- cluster.vcov(x.r, t.w56.r$GUNQID)
 x.mcse.r <- coeftest(x.r, x.clust.r)

 
 
  tr[[ref.names[i]]] <- x.mcse.r
  print(tr[[ref.names[i]]])
  
}


stargazer(tl[[ref.names[2]]], tr[[ref.names[2]]], tl[[ref.names[3]]], tr[[ref.names[3]]], tl[[ref.names[1]]], tr[[ref.names[1]]], tl[[ref.names[4]]], tr[[ref.names[4]]], covariate.labels = varlabels)



##Control Wave 6 and Control Wave 8 (Pre-Referendum: June 2016/ August 2017)

ones <- rep_len(1, length.out = nrow(control.sub))
twos <- rep_len(2, length.out = nrow(control.sub))


dvs <- Cs(closed.immigrants,hurt.standing.refugee, threaten.culture.refugee,overwhelm.welfare.refugee)

dvs.6 <- cbind(ones, control.sub$w8w6, control.sub$leaver, control.sub$closed.immigrants.6,
               control.sub$hurt.standing.refugee.6, control.sub$threaten.culture.refugee.6,
               control.sub$overwhelm.welfare.refugee.6,control.sub$GUNQID)

dvs.8 <- cbind(twos, control.sub$w8w6, control.sub$leaver, control.sub$closed.immigrants.8,
               control.sub$hurt.standing.refugee.8, control.sub$threaten.culture.refugee.8,
               control.sub$overwhelm.welfare.refugee.8,control.sub$GUNQID)


dvs.6 <- cbind(dvs.6)
dvs.8 <- cbind(dvs.8)


c.w68 <- as.data.frame(rbind(dvs.6, dvs.8))

names(c.w68) <- c("wave", "w8w6", "leaver",ref.names,"GUNQID")

c.w68.l <- subset(c.w68, leaver == 1)
c.w68.r <- subset(c.w68, leaver == 0)


cl <- vector(length(ref.names), mode = "list")
names(cl) <- ref.names

cr <- vector(length(ref.names), mode = "list")
names(cr) <- ref.names

#Leavers

for (i in 1:(length(ref.names))){
  
  modelformula <- paste(ref.names[i],right.side.eq)
  print(modelformula)
  
 x.l <- eval(substitute(lm(.modelformula, data = c.w68.l, 
                                           weights = w8w6), list(.modelformula = modelformula)))
  
 x.clust.l <- cluster.vcov(x.l, c.w68.l$GUNQID)
 x.mcse.l <- coeftest(x.l, x.clust.l)
 
 
 cl[[ref.names[i]]] <- x.mcse.l
 
 print(cl[[ref.names[i]]])
  
}

#Remainers
for (i in 1:(length(ref.names))){
  
  modelformula <- paste(ref.names[i],right.side.eq)
  print(modelformula)
  
  x.r <- eval(substitute(lm(.modelformula, data = c.w68.r, 
                                           weights = w8w6), list(.modelformula = modelformula)))
  
  
  x.clust.r <- cluster.vcov(x.r, c.w68.r$GUNQID)
  x.mcse.r <- coeftest(x.r, x.clust.r)
  
  cr[[ref.names[i]]] <- x.mcse.r
  
  print(cr[[ref.names[i]]])
  
}


stargazer(cl[[ref.names[2]]], cr[[ref.names[2]]], cl[[ref.names[3]]], cr[[ref.names[3]]], cl[[ref.names[1]]], cr[[ref.names[1]]], cl[[ref.names[4]]], cr[[ref.names[4]]], covariate.labels = varlabels)


